Table 2. Key model variables. 


Variable 

Formula 

Units 

Entrainment 

E = 2(L + W)av 

mV* 

Plume Water Volume Flux 

Q = LWv 

m J s" 

Momentum Flux 

M = LWv 2 

m 4 s* 

Temperature Flux 

O' 

II 

ft 

°C m 3 s' 1 

Dissolved Solids Flux 

F s =QSp w 

gs'* 

Dissolved O 2 and N 2 Fluxes 

F Di =QCi 

mol s' 

Gaseous O 2 and N 2 Fluxes 

Fqj =LWA?(v + v b )yi 

mol s 


Table 3. Non-linear differential flux equations. 


Water Volume Flux 

dz 

Momentum Flux 

dM _ Pa - p P gim 2 + p a - p w gLW (i_x, 2 ) 
dz p p Pp 

Temperature Flux 

- ex 

dz a 

Salinity Flux 

^£- = E Pa S a 
dz 

Dissolved Gas Flux 

dFoi =ECj + iH!!iK L (H j P i -Cj) 
dz v + v b 

Gas Flux 

dF Gi _ 4Tir 2 N KL(H _ p _ _ c .) 
dz v + v b 


Reservoir model 

CE-QUAL-W2 (W2) is a two-dimensional, laterally averaged, hydrodynamic and water quality model 
(USACE 1995). W2 includes horizontal, but not vertical, momentum, and accounts for momentum 
transfer from influent streams and the associated shear stress imparted to the surrounding water. W2 was 
applied to SHR for the year 1998 when the aeration system was operated. A temperature calibration was 
performed adjusting only the wind-sheltering coefficient. Based on the model calibration, no significant 
warming is expected to occur in the hypolimnion if the diffuser is not m operation. W2 was then set up to 
run during the aeration period. During diffuser operation there was no inflow into the reservoir and the 
only outflow was seepage loss and water supplied to the water treatment plant. For simplicity, all dissolved 
oxygen (DO) sinks were turned off, with DO treated as a conservative tracer. 


Coupled plume and reservoir model 

The location of the 360-m long diffuser changes in elevation with its deepest point being closest to the dam 
wall. Therefore, as shown in Figure 3, it is represented in W2 as 5 discrete segments m consecutive 
columns For the purpose of coupling, the plume model is used to compute the flow rate of ambient water 
entrained as the plume rises, and the flow rate, temperature, and oxygen concentration of the water that is 
detrained upon reaching the depth of maximum plume rise. In other words, entrainment removes water of 
known temperature and oxygen concentration from a range of depths, while detramment returns the sum of 
all the entrained water at the specified discharge location after adding the predicted amount of oxygen. As 
shown in Figure 3, this plume action is simulated within W2 by withdrawing water from the entrainment 



cells (marked with X in Figure 3) and discharging it at the “detrainment” cell (marked with 0 in Figure 3 ). 
Because of the changing temperature gradient during the course of the simulation, the plume rise height 
also changes. It was therefore necessary in some cases to specify a new detrainment cell during the course 
of the simulation, as indicated in Figure 3. 


Depth tmmts.) 



Although the reservoir is laterally averaged, the entrainment operation is adequately represented using a 
simple withdrawal. However, the entire plume discharge is mixed into the detrainment cell. Because the 
theoretical plume width at the point of detrainment approaches infinity, this instantaneous mixing appears 
reasonable to a first approximation. However, the selection of detrainment cell size might be an important 
consideration when the coupling procedure is refined. 


The “manual” coupling procedure was handled as follows. First, the initial reservoir conditions on Sep. 28 
were used in the plume model to obtain the flow rate of entrained and detrained water This^procedure was 
performed for each of the 5 diffuser segments shown in Figure 3. Then, starting on Sep. 28 at 12h00, and 
using the 5 sets of predicted data, the coupled plume/W2 model was run for 6 hours. The predicted 
ambfent temperature and DO profiles obtained at the end of this 6-hour period were then used asmputto 
the plume model to generate a new set of predicted plume entrainment and detrainment data for each of the 
5 plume segments. The entire procedure was repeated until Oct. 9 at 12h00. 


RESULTS AND DISCUSSION 

The results of the preliminary coupling procedure are shown in Figures 4 and 5. Excluding the water below 
the diffuser. Figure 4 shows that the coupled model predicts mixmg and warming induced by plume 
operation quite accurately. The thermocline is eroded during the entire period of operation even though 
te'C’olSy partially^penetrates the thermocline during the first day of 

diffuser segments shown in Figure 3. Thereafter, the plume stops between 1 and 3 meters below the 
thermocline The results suggest that the flow of detrained water causes sufficient mixing nea t 
thermocline to result in almost exactly the right amount of wanning in the hypolimmon. Factors causing 










mining include vertical dispersion of horizontal momentum, momentum transfer from the detraining plume 
water, and horizontally generated shear. A more thorough investigation into the dominant mixing 
mechanisms is required. 

Although all oxygen sinks are turned off in the W2 model, Figure 5 shows that the coupled model predicts 
the evolution of hypolimnetic DO fairly well. Once the effects of sediment and other oxygen demands are 
correctly incorporated in W2, the predicted oxygen concentrations will not increase as rapidly. However, 
as currently represented, the plume model under-predicts the rate of oxygen addition because oxygen 
transfer ceases when the plume detrains. Bubbles are still present in the hypolimnion and will continue to 
rise through the water column. In effect, this produces a weaker secondary plume that results in some 
additional mixing and oxygen transfer in the hypolimnion, as well as at higher locations in the reservoir. 
The effects of the secondary plume will be investigated, and if significant relative to the primary plume, 
will be incorporated in the coupled model. 



Temperature (°C) 

Figure 4. Observed versus predicted 
hypolimnetic temperature. 


Dissolved Oxygen (g rrf 3 ) 

Figure 5. Observed versus predicted 
hypolimnetic DO. 


FUTURE WORK 

While the preliminary results are very encouraging, the iterative “manual” coupling procedure is time 
consuming. The plume model will therefore be incorporated as a subroutine in W2 so that the computation 
can be done automatically. Further testing of the coupled model will be carried out on new data collected 
at SHR in Virginia as well as at Lake Baldegg in Switzerland. Ultimately, the coupled model could also be 
used to improve the calibration of the plume model. 
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[i] An existing linear bubble plume model was improved, and data collected from a 
full-scale diffuser installed in Spring Hollow Reservoir, Virginia, were used to validate 
the model. The depth of maximum plume rise was simulated well for two of the three 
diffuser tests. Temperature predictions deviated from measured profiles near the maximum 
plume rise height, but predicted dissolved oxygen profiles compared very well with 
observations. A sensitivity analysis was performed. The gas flow rate had die greatest effect 
on predicted plume rise height and induced water flow rate, both of which were directly 
proportional to gas flow rate. Oxygen transfer within the hypolimnion was independent 
of all parameters except initial bubble radius and was inversely proportional for radii 
greater than approximately 1 mm. The results of this work suggest that plume dynamics 
and oxygen transfer can successfully be predicted for linear bubble plumes using the 
discrete-bubble approach. 

Citation: Singleton, V. L„ P. Gantzer, and J. C. Little (2007), Linear bubble plume model for hypolimnetic oxygenation: Full-scale 
validation and sensitivity analysis. Water Resour. Res., 43, W02405, doi:10.1029/2005WR004836. 


1. Introduction 

[ 2 ] Bubble plumes are used in a variety of industrial and 
environmental applications including mixing in chemical 
reactors, stripping of dissolved gases, containment of spills, 
prevention of ice formation, protection of harbors from 
damaging waves \Fannel 0 p et al, 1991], and destratifica¬ 
tion of lakes and reservoirs [Schladow, 1992]. In addition to 
airlift aerators [Burris et al., 2002] and Speece cones 
[McGinnis and Little, 1998], bubble plumes are commonly 
used for hypolimnetic aeration and oxygenation, which 
preserves stratification of water bodies while adding oxygen 
to the deepest layer. Hypolimnetic anoxia negatively affects 
the drinking water treatment process, cold-water fisheries, 
and water quality downstream of hydropower reservoirs. In 
the United States, releases from hydropower reservoirs 
typically must comply with state water quality criteria for 
minimum dissolved oxygen (DO) concentrations [Peterson 
et al, 2003], Oxygen depletion may lead to increases in 
hydrogen sulfide, ammonia, and phosphorus and the release 
of reduced iron and manganese from the sediments. Hydro¬ 
gen sulfide, iron, and manganese in drinking water usually 
require additional treatment [Cooke and Carlson, 1989], 
Finally, hypoxia can affect sex differentiation and develop¬ 
ment in fish, resulting in male-dominated populations with 
reduced reproductive success [Shang et al., 2006]. 

[ 3 ] A bubble plume model to predict oxygen transfer 
from linear diffuser systems was presented by McGinnis et 
al. [2001], based on the model for a circular diffuser 
developed earlier by Wiiest et al. [1992]. While several 
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models for point-source or circular bubble plumes have 
been proposed [Asaeda and Imberger, 1993; Brevik and 
Kluge, 1999; Ditmars and Cederwall, 1974; Fannelep and 
Sjoen, 1980; Johansen, 2000; Kobus, 1968; McDougall, 
1978; Milgram, 1983; Rayyan and Speece, \911\Sahoo and 
Luketina, 2003; Schladow, 1992; Wiiest et al., 1992; Zheng 
et al., 2002], less work has been conducted on linear (also 
referred to as line, two-dimensional, or planar) bubble 
plumes. Kobus [1968] developed one of the first analytical 
models for linear bubble plumes, which uses an empirical 
correlation to calculate buoyancy flux. Ditmars and Cederwall 
[1974] presented a model similar to that of Kobus [1968] but 
included bubble slip velocity. Brevik [1977] proposed a 
phenomenological theory for two-dimensional bubble 
plumes comparable with that of Ditmars and Cederwall 
[1974], except that kinetic energy was used to predict 
entrainment. Wilkinson [1979] proposed that full-scale 
linear plumes could be characterized by a Weber number. 
Laureshen and Rowe [1987] presented a model for two- 
dimensional bubble plumes in which plume spreading, 
entrainment, and momentum amplification were assumed 
to be functions of the plume Weber number and empirical 
constants. Fannelep et al. [1991] developed a model for 
linear plumes in shallow water and studied the resulting 
surface currents and recirculation cells. Last, Brevik and 
Kluge [1999] expanded an existing model for linear bubble 
plumes to account for vertical turbulence. Although much 
insight into plume dynamics was gained, none of these 
models for linear or two-dimensional bubble plumes 
accounted for ambient stratification or gas transfer. The 
first linear bubble plume model to include gas transfer was 
presented by McGinnis et al. [2001], who converted the 
circular bubble plume model of Wiiest et al. [1992] to linear 
geometry. The incorporation of gas transfer is critical 
because the rapid dissolution rate of oxygen, and nitrogen 
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when compressed air is used, strongly influences the buoy¬ 
ancy of the plume [Wiiest et al., 1992]. Gas transfer is 
especially important in deep water bodies and for weak 
plumes because the increased contact time allows greater 
gas exchange. Last, the prediction of oxygen addition from 
hypolimnetic oxygenation systems is facilitated. Despite 
the usefulness of the linear bubble plume model, it has not 
yet been validated at full scale and over a range of operating 
conditions. 

[ 4 ] Using extensive, high-spatial-resolution conductivity 
and temperature as a function of depth (CTD) transect data 
collected in Spring Hollow Reservoir (SHR), Virginia, 
during diffuser operation in 2003 and 2004, the perfor¬ 
mance of the linear bubble plume model is evaluated. The 
motivation for this work includes verification of model 
performance prior to use for design and investigation of 
critical model parameters through sensitivity analysis. Also, 
the accuracy of model predictions for depth of maximum 
plume rise (DMPR) and induced water flow rate should be 
assessed prior to coupling with lake/reservoir hydrody¬ 
namic and water quality models, such as CE-QUAL-W2 
[McGinnis et al., 2001]. In this paper an improved linear 
bubble plume model is presented, observations and model 
predictions are compared, and results of a sensitivity 
analysis are discussed. 

2. Bubble Plumes in Stratified Water Bodies 

[ 5 ] During hypolimnetic oxygenation with bubble 
plumes, compressed gas is continually supplied to diffusers, 
usually located immediately above the sediments, and is 
allowed to bubble freely. A gas-water plume mixture that is 
less dense than the ambient water is created, which causes 
the mixture to ascend through the hypolimnion. As the 
mixture rises, ambient water is entrained into the plume, and 
the plume width increases. The entrained fluid produces a 
double-plume structure, consisting of an inner core that 
contains the bubble-water mixture surrounded by an outer 
annulus that contains plume water relatively free of bubbles 
[McDougall, 1978]. As the outer annulus entrains stratified 
hypolimnetic water, the plume width increases and the 
density decreases. When the negative buoyancy of the 
entrained fluid exceeds the positive buoyancy imparted by 
the bubbles, the plume detrains water at a rate nearly equal 
to that previously entrained [ Lemckert and Jmberger, 1993]. 
At this depth, the velocity of the relatively dense water 
within the plume decreases to zero, and the plume stops 
rising. The detraining plume water then forms an annular 
downward flow immediately outside the outer annulus of 
the upward flowing plume water. The detraining plume 
water entrains ambient water until a depth of neutral buoy¬ 
ancy is reached, where a horizontal intrusion is created into 
the hypolimnion [Asaeda and Imberger, 1993], The undis¬ 
solved bubbles remaining in the bubble-water mixture 
separate from the inner core flow and continue to rise to 
the surface, repeating the entire process. 

3. Linear Bubble Plume Model 

[6] The linear bubble plume model utilizes the discrete- 
bubble approach, which has also been applied to the airlift 
aerator and Speece cone and was recently reviewed in detail 
by Singleton and Little [2006]. The linear bubble plume 


model is composed of horizontally integrated equations 
based on the conservation of mass, momentum, and heat. 
Eight flux equations are solved simultaneously to predict 
water flow rate, plume temperature, oxygen and nitrogen 
transfer and concentration, salinity, and plume rise height, 
given diffuser geometry and depth, applied gas flow rate, 
and initial bubble size (Tables 1 and 2). The model accounts 
for density stratification due to vertical temperature and 
salinity gradients. Entrainment is assumed to be propor¬ 
tional to the local (with respect to depth) plume water 
velocity and perimeter. Bubble size varies as the bubbles 
rise due to expansion and dissolution, and bubble slip 
velocity and gas transfer coefficients are functions of 
bubble radius [Wiiest et al, 1992], Also, Henry’s constants 
for oxygen and nitrogen are functions of temperature 
[Wiiest et al., 1992], The bubble plume model equations 
were originally developed by Wiiest et al. [1992] for 
circular geometry but were modified by McGinnis et al. 
[2001] for the linear geometry of the system installed in 
SHR. The equations that include the spreading coefficient 
(A) were recently refined [Singleton and Little, 2005] to 
more accurately reflect the geometry of the plume at the 
ends of the linear diffuser (Tables 1 and 2), which is 
approximated in plan view as a long, thin rectangle. 
Additional refinements to the model of McGinnis et al. 
[2001], which are detailed in the following paragraphs, 
include use of a correlation to calculate initial bubble size 
[McGinnis and Little, 2002], correction of the entrainment 
coefficient (a) and A for top-hat profiles, and use of a 
Froude number (Fr) to calculate initial water velocity 
[Fischer et al., 1979; Wiiest et al., 1992]. Also, water 
quality profiles from the plume near-field, as opposed to 
the reservoir far-field, were used as boundary conditions 
[McGinnis et al., 2004]. 

[7] The entrainment coefficient and X were set at 0.11 and 
0.93, respectively. These values were derived by Fannelop 
et al. [1991] by fitting Gaussian profiles to laboratory data 
and were modified for the top-hat profile assumption of the 
model using [Fannelop and Sjsen, 1980] 

ct r = \fiac (1) 

= (2) 

where the subscripts T and G refer to top-hat and Gaussian 
profiles, respectively. Ffrr simplicity, top-hat or uniform 
profiles are assumed for water velocity, temperature, 
salinity, dissolved and gaseous constituents, and bubble 
velocity [Wiiest et al., 1992]. Other model assumptions are 
as follows: (1) The linear plume width W for temperature 
and dissolved constituents is equal to the width of the plume 
velocity profile, whereas the bubbles are confined to an 
inner core of width A W (A < 1); (2) ambient currents are 
negligible; (3) the diffuser produces bubbles at a constant 
rate and uniform size that are evenly distributed over the 
cross section of initial width A W 0 ; (4) bubble coalescence is 
neglected; (5) initial water properties of the plume are those 
of ambient water at the diffuser depth; and (6) exchange of 
gases other than oxygen and nitrogen is not considered. 

[8] The model predictions are strongly dependent on the 
initial plume conditions and the plume boundary conditions 
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Table 1. Key Variables of the Linear Bubble Plume Model 2 


Variable Formula Units 


Entrainment factor E = 2(L + W)av mjs 

Plume water volume flux Q = LWv m /s 

Momentum flux M - LWv 2 ^ 

Temperature flux Ft — QJ P °C m /s 

Dissolved solids flux F, = QSp„ Itgfs 

Dissolved 0 2 and N 2 fluxes Fd, = QC, mol/s 


Gaseous 0 2 and N 2 fluxes F a , - - W (1 - A)] (v + v b )y, mol/s 

“Revised from McGinnis el al. [2001]. 


of temperature, dissolved oxygen, and salinity. Initial con¬ 
ditions for the bubble plume model were determined as 
detailed by Wiiest et al. [1992], except for the following 
deviations. Initial bubble size is calculated using the corre¬ 
lation developed by McGinnis and Little [2002] for the type 
of linear diffiiser installed in SHR: 

d X2 = 1.12 + 0.938?. (3) 

The correlation was determined using measured Sauter- 
mean bubble diameter ( d 3 , 2 ) values of 1.1-2.2 mm 
collected over actual unit gas flow rates q of 0.08- 
0.88 m 2 /h at the diffuser. 

[ 9 ] For the circular bubble plume model, Wiiest et al. 
[1992] proposed that the induced vertical water velocity v 
at the diffuser depth is equivalent to the initial plume 
water velocity. To estimate the initial velocity, Wiiest et al. 
[1992] defined a densimetric Fr and utilized a relationship 
between the local Richardson number (Ri) and Fr derived 
by Fischer et al. [1979] for single-phase, round buoyant 
jets discharging vertically. The circular bubble plume 
model was recently validated by McGinnis et al. [2004], 
so a similar procedure was used to determine the initial plume 
velocity for the linear bubble plume. Corresponding equa¬ 
tions for Fr were derived for planar or linear plumes using 
relationships presented by Fischer et al. [1979] to obtain 

Fr = Ri~‘ (4) 


the source where the flow is more like a plume [Fischer et 
al., 1979]. Consequently, Fr for a planar or linear plume is 
equal to 1.26, except close to the source. Wiiest et al. [1992] 
assumed that the bubble slip velocity near the source was 
relatively low, so the initial Froude number [Fr 0 ) should be 
equal to the value for a single-phase plume. Unlike the 
circular plume model, the Fr profiles predicted by the linear 
bubble plume model continually increase with depth for the 
diffiiser installed in SHR (not shown). For Fr 0 < 2.0, the 
plume velocity initially increases with decreasing depth 
immediately above the linear diffiiser (Figure la). This 
effect was also predicted by Farmelep and Webber [2003] 
for buoyant plumes rising from areal sources, where a point 
of maximum velocity occurred above the source. Addition¬ 
ally, the plume neck (point of minimum radius or width) is 
always below the point of maximum velocity in the plume 
[Fannelep and Webber, 2003]. A neck is not predicted for 
the linear bubble plume when Fr a = 2.0 (Figure lb). 
However, necking or contraction will only occur when the 
momentum immediately above the source is relatively low 
and/or the entrainment coefficient is relatively low 
[Fannehsp and Webber , 2003]. Also, Wiiest et al. [1992] 
reasoned that a plume from an open source (diffiiser above 
sediments) may not contract because initially entrained 
water is not obstructed. Therefore Fr a for an open source 
will likely be higher than that for a closed source (diffuser 
resting on sediments). Because the linear diffuser in SHR 
is an open source (Figure 2), Fr„ for the linear bubble 
plume was assumed to be 2.0, and a sensitivity analysis 
was conducted to determine the effect of varying Fr„. 

[ 10 ] The differential flux equations of the linear bubble 
plume model (Table 2) were solved numerically using the 
fourth-order Runge-Kutta method. Further information on 
the general solution procedure, equations of state, and 
model assumptions is provided by Wiiest et al. [1992] and 
McGinnis et al. [2004]. The model calculations are only 
valid over the plume rise height, up to the DMPR. When the 
plume stops rising, a secondary plume may form above as 
bubbles that are not completely dissolved continue to rise 
[Asaeda and Imberger, 1993; McDougall, 1978; Schladow, 
1992], This phenomenon can occur when a bubble plume is 
released into strong density stratification. 


rr = -;-J , w 

\>~Wg{p a - P p )/Pp]~ 

where g is gravitational acceleration, p„ is the ambient water 
density, and p p is the bubble plume density. The local Ri for 
planar jets has a constant value of 0.735 at distances from 


4. Application to Spring Hollow Reservoir, 
Virginia 

4.1. Field Data Collection 

[ 11 ] To fully evaluate the linear bubble plume model, 
experimental data for boundary conditions, rise height, and 
in-plume constituent profiles are required. Testing was 


Table 2. Nonlinear Differential Flux Equations of the Linear Bubble Plume Model" 


Flux 

Equation 

Water volume 

Momentum 

Temperature 

Salinity 

Dissolved gas 

Gas 

dQldz = E 

dMIdz = [(p„ - p„yp p lgEW + l(p» - PpyPplgM v - - *)J 

rfFj/rfe = ET a 
dFJdz = E Pa S a 

dF n jdz = EC„ + [AizFNK.v + v h )]KdH,P, - C,) 
dFgjdz = — [4m- 2 W/(v + v*)] K L (fijPi - C,) 


"Revised from McGinnis et al [2001]. 
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Figure 1. Effect of initial Froude number (Fr„) on 
predictions of plume water velocity and plume width using 
linear bubble plume model. 


conducted using a full-scale linear diffuser (Figure 2) 
installed in Spring Hollow Reservoir, Virginia (Figure 3). 
Constructed in 1995, SHR is a small monomictic, mesotro- 
phic side-stream reservoir that is generally stratified from 
May to December. The reservoir is managed by the Western 
Virginia Water Authority and serves as one of the principle 
drinking water sources for Roanoke County. The water 
body has a maximum depth of 65 m and a maximum 
surface elevation of 431 m. The approximate surface area 
and volume are 0.54 km 2 and 12.4 x 10 6 m 3 , respectively. 
To prevent anoxia in the hypolimnion and the associated 
deterioration of raw water quality, a linear diffuser equipped 
with fine-bubble porous hoses was installed in 1997 
(Figure 3). The 305-m-long diffuser can be supplied with 
compressed air or pure oxygen at various gas flow rates and 
is located in the deepest portion of the reservoir (368- 
372 m elevation). On the basis of an average surface 
elevation of 430 m, the depth of the diffuser during testing 
ranged from 58 to 62 m along its length. Table 3 provides 
the diffuser operating parameters in 2003 and 2004. 

[ 12 ] Diffuser tests were performed in 2003 using com¬ 
pressed air (21% 0 2 ) supplied at a high gas flow rate 
(45 Nm 3 /h average) during 29 June to 14 July and pure oxygen 
(97% 0 2 ) supplied at a low flow rate (II Nm 3 /h average) 
during 14-26 August. A third test was conducted in 
2004 using pure oxygen, but at a higher gas flow rate 
(40 Nm 3 /h maximum) during 22 October to 5 November. In 
2003, the diffuser quickly mixed the rather small hypolim- 
netic volume during both tests, and the water quality con¬ 
ditions on the dates of data collection (2 July and 17 August 
for compressed air and pure oxygen, respectively) were by 
that time relatively homogeneous as a result of plume- 
induced mixing. One of the primary objectives of the 2004 


experiments was to maximize the plume signature in the 
hypolimnion and to increase confidence in the linear 
plume model validation under a different set of boundary 
conditions. In 2004, the data were therefore collected on 
24 October, soon after start of diffuser operation. Addition¬ 
ally, the 2004 test was performed later in the stratified season 
to maximize the ambient DO and temperature gradients in 
the hypolimnion. 

[ 13 ] To establish appropriate boundary conditions for the 
plume model, characterization of the plume near-field 
environment is necessary [McGinnis et al., 2004]. Therefore 
the data collected included numerous high-spatial-resolution 
CTD (Sea-Bird model SBE 19plus; 4-Hz sampling rate) 
transects measured almost daily before, during, and after 
diffuser operation. The CTD profiler was also equipped 
with a DO probe (1.4-s response time measured at 20°C). 
Profiles were obtained laterally across the diffuser at 0.5-m 
increments for 0-10 m, 2-m increments for 10-20 m, and 
5-m increments for 20-40 m from the centerline of the 
diffuser in both directions (Figures 3 and 4). (The diffuser 
centerline location is shifted to the left for 2004 (Figures 4e 
and 4f) because the diffuser was repositioned earlier in the 
year. Also, the operational length of the diffuser was 
decreased for that year (Table 3).) 

[u] The geometry of SHR affects the extent and rate of 
circulation within the hypolimnion induced by the bubble 
plume. Because of the relatively small size of SHR, oper¬ 
ation of the linear diffuser created uniform conditions below 
the thermocline within days after startup. While SHR 
bathymetry influences plume-induced mixing in the hypo¬ 
limnion, the effect on short-term plume operation is negli¬ 
gible because the time that individual bubbles spend in the 
hypolimnion is of the order of minutes. The bubbles and 
resulting plume experience a pseudo steady state with 
respect to ambient conditions. Effects on data collection 
due to the ends of the linear diffuser were assumed to be 
negligible because the lateral profile location was over 
150 m from a diffuser end (Figure 3), and the diffuser is 
designed to release a uniform gas flow along its length. 

4.2. Observations and Model Validation 

[ 15 ] Plume rise height, spreading, and constituent profiles 
predicted by the linear bubble plume model were compared 
with experimental observations. Critical model input param¬ 
eters for the three test conditions (2 July 2003, 17 August 
2003, and 23 October 2004) are shown in Table 3. The 
boundary conditions were obtained from averaged near¬ 
field lateral profiles [McGinnis et al., 2004] (±2 m and ±1 m 
from plume centerline for 2003 and 2004, respectively) and 
differed significantly between 2003 and 2004 (Figure 5). 
(The presence of two thermoclines in 2003 is due to the 
pumped reservoir inflow that discharges at 396-m elevation, 
which corresponded to approximately 34-m depth during 
diffuser testing. The lower thermocline delineates the effec¬ 
tive hypolimnion for the oxygenation system.) 

[ 16 ] Measured contours of temperature and DO are 
shown in Figure 4, along with corresponding model pre¬ 
dictions for plume width and the DMPR. The actual plume 
boundaries are not well defined in the contour plots, so 
comparison with predicted plume widths is difficult. The 
lack of distinct plume boundaries was due to the almost well 
mixed conditions in the hypolimnion as a result of diffuser 
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Figure 2. Photograph and schematic of linear bubble 
plume diffuser in Spring Hollow Reservoir, Virginia. 
Courtesy of Mark Mobley, Mobley Engineering, Inc. 


operation, particularly for 2003 (Figures 4a—4d). Also, use 
of compressed air did not produce a strong DO plume 
signature for the July test compared with the August and 
October tests with pure oxygen. The actual plume rise 
height is easier to distinguish, especially in the DO contours 
for August and October. The predicted depths of maximum 
plume rise are 38.2, 46.1, and 48.8 m for July, August, and 
October, respectively. For July and October, the DMPR is 
simulated well by the model (Figures 4a, 4b, 4e, and 4f). 
However, the model appears to underestimate the plume rise 
height for August, when the gas flow rate was compara¬ 
tively low (Table 3). The under-predicted DMPR may have 
been due to an overestimated value for a. In a detailed study 
of round plumes, Milgram [1983] found that a is directly 
proportional to the plume gas holdup or fraction. However, 
the linear bubble-plume model assumes that a is constant 
(Table 3). 

[n] The structure of the plumes is similar to those 
observed by Asaeda and Imberger [1993] for round bubble 
plumes in weak stratification (Figure 4). Depending on the 
gas flow rate and stratification strength, three types of 
horizontal intrusions from the plume were reported. The 
pattern of the DO contours for October (Figure 4f) closely 
resembles type 1, which corresponds to a high gas flow rate 
or weak stratification [Asaeda and Imberger , 1993]. This is 
similar to a single plume impinging on a free surface, which 


is analogous to the thermocline in SHR. Although the 
plume for July is not easily discerned from the temperature 
or DO contours, the structure is most likely similar to a type 
1 plume. The July plume appears to detrain primarily at the 
lower thermocline with one strong intrusion, as evidenced 
by accumulation of higher oxygenated water near the top of 
the plume (Figure 4b). The plume for August is best 
classified as type 3, which is for low gas flow rates or 
strong density stratification [ Asaeda and- Imberger , 1993]. 
Type 3 plumes do not have steady intrusions, but instead are 
characterized by alternating, collapsing eddies that cause the 
plume to meander. 

[is] Referring to the October test (Figure 4f), the higher 
DO concentrations at lower depths adjacent to both sides of 
the plume were likely the result of detrained water that sinks 
past the equilibrium depth due to momentum. The temper¬ 
ature isotherms were also depressed immediately beside the 
plume (Figure 4e). This phenomenon was also observed by 
McGinnis et al. [2004]. The DO concentration immediately 
above the predicted DMPR and near the vertical plume 
centerline for October 2004 is higher than the ambient 
concentration (Figure 4f). This could have been earned by 
the formation of a secondary plume above the DMPR 
resulting from incompletely dissolved bubbles. The model 
estimates that the bubble size at the top of the first plume 
was about 5 x 10" 5 m for October 2004, which is relatively 
large compared with July 2003 and August 2003 (Figure 6f). 
These undissolved bubbles could have created a secondary 
plume, which entrained oxygenated water from the detrain- 
ment of the first plume and carried it higher into the water 
column. 

[ 19 ] Vertical profiles of constituents and properties within 
the plumes were also predicted for the three diffuser tests 
(Figure 6). For the July test with compressed air, the higher 
gas flux creates a greater buoyancy flux and a higher initial 
plume water velocity (Figure 6d). Also, the concentration 
driving force for oxygen transfer is lower compared with 
pure oxygen, which decreases the rate of bubble dissolution 
with depth (Figure 6f). These effects result in a higher 
DMPR for July compared with August (Figures 4 and 6). 
The model predictions for 2004 differ from those of 2003 
because of the differing boundary conditions (Figure 5). The 
plume rise height for the October test with pure oxygen is 
less than that for August, even though the gas flux was more 
than tripled (Table 3). The ambient temperature, and hence 
density, stratification was stronger in 2004, which provided 
greater negative buoyancy to decrease plume momentum. 
The stronger ambient density stratification also caused the 
plume velocity to decrease more rapidly with depth in 
2004 despite a higher initial velocity from the diffuser 
(Figure 6d). The lower plume rise height in October resulted 
in lower plume water flow rates than July (Figure 6e), even 
with comparable gas fluxes applied (Table 3). 

[ 20 ] Average temperature, DO, and density profiles 
within the plumes were also measured (Figures 6a-6c). 
The temperature, and consequently density, predictions for 
July and October deviate from the measured profiles where 
the plumes reach the top of the hypoiimnion (Figures 6a 
and 6c), or where the rate of plume spreading is greatest 
(Figure 4). The model underpredicts the final plume 
temperature by approximately 0.3° and 0.2°C for July 
and October, respectively. One reason for the discrepancy 
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may be that the mode! assumes a is constant In addition 
to the dependence on plume gas holdup, Milgram [1983] 
also found that a for a circular plume is directly propor¬ 
tional to the local plume radius. As the plume approaches 
its maximum rise height, the width increases rapidly 


(Figure 4). If linear plume dynamics are similar to circular 
plumes, then a for linear plumes may also increase as the 
plume width increases. Additionally, the boundary profiles 
selected may not accurately reflect the actual ambient 
conditions immediately adjacent to the plume along its 
entire rise height Measured vertical CTD profiles at the 
estimated plume width were averaged and used for boundary 
conditions. However, the plume width varies greatly with 
depth (Figure 4), so use of vertical profiles at a single lateral 
distance from the diffuser is not appropriate. Also, the 
average plume width was visually estimated from the 
temperature and DO contour plots (Figure 4), but the actual 
plume boundaries are not well defined. 

[ 21 ] The model predicts the plume DO profiles well for 
all three diffuser tests (Figure 6b). For July and August the 
model characterizes the initial increase in DO immediately 
above the diffuser quite accurately. The initial rapid increase 
in DO for August and October is due to the higher oxygen 
saturation concentration at depth with the use of pure 
oxygen. The high hydrostatic pressure causes the oxygen 
transfer rate to be almost independent of the ambient DO 
concentration at the depth of diffuser. By contrast, the initial 
DO increase for the July test with compressed air is more 
modest (Figure 6b). The shape of the predicted DO profile 
for October differs somewhat from the experimental data, 
even though the final values at the DMPR differ by only 
0.3 g/m 3 . The hypolimnion and hence the plume near-field 
were more heterogeneous in 2004 than 2003, which may 
have contributed to the overprediction at lower depths if 
the selected boundary profiles did not accurately represent 
water entrained into the plume. 

[ 22 ] Another source of inaccuracy could be the correla¬ 
tion equation used to calculate initial bubble size (equation 
(3)). This relationship was developed using data collected 
over actual air flow rates per unit length of diffuser of 0.08- 
0.88 m 2 /h [McGinnis and Little , 2002], The actual gas flow 
rates per unit length of diffuser for the July, August, and 
October tests were 0.019, 0.0063, and 0.024 m 2 /h, respec¬ 
tively. Therefore the initial bubble diameters used in the 
model were extrapolated beyond the valid correlation range. 


Table 3. Conditions for Linear Bubble Plume Model Validation and Sensitivity Analysis for Spring Hollow Reservoir, Virginia 2 


Parameter 

2 July 2003 

17 Aug 2003 

23 Oct 2004 

Sensitivity Range 

Oxygen in gas supply, % 

21 

97 

97 

n/a 

Entrainment coefficient [ ] 

0.11 

0.11 

0.11 

0.05-0.2 

Spreading coefficient [ ] 

0.93 

0,93 

0.93 

0.5-1.0 

Initial Froude number [ ] 

2.0 

2.0 

2.0 

1.0-3.0 

Operational diffuser length, m 

300 

300 

250 

60-625 

Initial plume width, m 

0.16 

0.16 

0.16 

n/a 

Initial plume area, m 2 

50 

50 

42 

10-100 

Gas flow rate, b Nm 3 /h 

38 

13 

40 

1-400 

Initial gas flux,' m/h 

0.76 

0.26 

0.96 

0.02-9.5 

Initial bubble radius, m 

5.7 x 10~ 4 

5.6 x 10 -4 

5.7 x 10 4 

10 -4 — J0“* 

Diffuser depth, d m 

60 

60 

59 

n/a 

Reservoir maximum depth, m 

64 

65 

64 

n/a 

Reservoir surface area, 10" m 2 

0.53 

0.53 

0.53 

n/a 

Reservoir total volume, 10 6 m 3 

12 

12 

12 

n/a 


“Reservoir conditions on testing days are also included. 
b Nm 3 denotes I m 3 of gas at 1 bar and 0°C. 

'Sensitivity analysis range refers to varying gas flow rate while maintaining constant initial plume area. 
‘'Depth at location of lateral CTD transect. 
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Figure 4. Measured (left) plume temperature (°C) and (right) dissolved oxygen (DO) (g/m 3 ) contours 
with linear bubble plume model predictions for diffuser operation with air (2 July 2003) and pure oxygen 
(17 August 2003 and 23 October 2004) in Spring Hollow Reservoir, Virginia. Contours were interpolated 
from CTD profiles collected at locations indicated by small black squares along the bottom of each plot 


4.3. Sensitivity Analysis 

[ 23 ] A sensitivity analysis was performed with the linear 
bubble plume model to determine the effects on plume rise 
height, oxygen transfer efficiency, and induced water flow 
rate, because these parameters are important for design and 
operation of hypolimnetic aeration/oxygenation systems. 
Parameter perturbation was used, in which input variables 
are independently adjusted to determine their individual 
effects on model predictions. The model variables investi¬ 
gated are either difficult to measure or can be controlled 
through system design or operation (Table 3). Currently, a 
and A for the linear plume model are empirical constants 
\Fannel 0 p et al, 1991], and the initial plume water velocity 
is calculated using a densimetric Fr [Wiiest et al, 1992]. 
The initial plume area, gas flow rate, and, to a lesser extent, 
initial bubble size can be controlled through diffuser design 
and operation. The sensitivity of model predictions to 


ambient dissolved nitrogen was also examined for the 
standard case using compressed air (2 July 2003). The 
model assumes that the background dissolved nitrogen 
concentration is equivalent to the saturated value at atmo¬ 
spheric partial pressure and the average hypolimnetic water 
temperature. 

[ 24 ] The DMPR is most influenced by gas flow rate and 
initial bubble radius (Figures 7e and 7f). The gas flow rate is 
directly related to the density of the plume bubble-water 
mixture and subsequently the positive buoyancy and up¬ 
ward momentum of the plume (Table 2). For gas flow rates 
greater than approximately 100 Nm 3 /h, the plume rise 
height is controlled more by the depth of the thermocline 
in SHR (Figure 5a). As the plume rise height approaches the 
thermocline, further increases in the buoyancy and momen¬ 
tum fluxes can not overcome the strong ambient density 
stratification. The 2003 predictions are more a function of 
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Figure 5. Input boundary conditions for linear bubble plume model validation and sensitivity analysis. 
Data collected from Spring Hollow Reservoir, Virginia, during diffuser operation with compressed air 
(2 July 2003, crosses) and pure oxygen (17 August 2003, solid circles, and 23 October 2004, open circles). 


gas flow rate than those for October 2004 because of 
differing boundary conditions (Figure 5). 

[ 25 ] The plume rise height is moderately sensitive to 
initial bubble radius (Figure If). For initial radii less than 
about 1 mm, the DMPRs are nearly independent of bubble 
size. As initial bubble sizes increase, the plumes ascend 
higher and reach a maximum height for a radius of about 
6 mm for each diffuser test. The shapes of the curves can 
be attributed to the dependence of bubble rise velocity 
and gas transfer coefficients on bubble radius [Wiiest et 
al., 1992], The DMPR predictions for August are more 
sensitive to initial bubble Tadius than those for July and 
October (Figure 7f). The plume on 17 August 2003 did 
not have sufficient buoyancy and momentum to rise to the 
thermocline because a relatively low gas flow rate was 
applied (Figure 4c and Table 3). By contrast, the plume 
dynamics for July and October were influenced to a greater 
degree by the thermocline, which in effect damped the 
sensitivity of plume rise height to the initial bubble radius. 

[ 2 s] The predicted DMPR for the linear bubble plume is 
virtually independent of Fr a and A over the ranges analyzed 
(Figures 7b and 7c). The plume rise height is moderately 
sensitive to ct Entrainment into the plume is a function of 
plume size, water velocity, and a (Table 1), and entrainment 
of ambient water decelerates the plume. Plume rise is 
influenced to a somewhat greater degree by the initial plume 
area. As the plume area is increased, the buoyancy flux 
decreases because the gas flow rate is constant. Plume rise 
height was unaffected by variations in the ambient dissolved 
nitrogen concentration from 50 to 200% saturation in the 
hypolimnion for the 2 July 2003 diffuser test. Similar to 
the analysis for DMPR, induced water flow rate at the top 
of the plume was found to be insensitive to Fr 0 and A and 


somewhat more sensitive to a and the initial plume area 
(results not shown). Plume water flow rate is most 
influenced by gas flow rate and initial bubble radius. 
Higher gas flow rates produce greater buoyancy and 
momentum fluxes, which results in greater plume rise 
heights ahd increased entrainment. 

[ 27 ] Oxygen transfer efficiency (total mass of oxygen 
transferred relative to initial mass of oxygen in bubbles) 
within the hypolimnion was independent of all parameters 
except initial bubble radius, decreasing from nearly 100% to 
around 20% for air and pure oxygen as the initial radii 
increased from approximately 1 mm to 1 cm (Figure 8). 
Oxygen transfer can continue above the DMPR if the 
bubbles are not dissolved. However, secondary plumes are 
not accounted for in the model. Even though undissolved 
bubbles at the top of the plume may continue to transfer 
oxygen during ascent, the oxygen may not be added at the 
desired depth (i.e., below the thermocline). A local mini¬ 
mum with respect to induced water flow rate at the top of 
the plume is predicted for a bubble radius of about 1 mm for 
the conditions in SHR (not shown). This suggests that while 
maximum oxygen transfer efficiency can be achieved with a 
1-mm initial bubble radius, vertical water circulation, and 
hence oxygen distribution in the hypolimnion, will not be 
optimized. As the bubble radius is increased from 1 mm, 
induced water flow rate increases as well for radii up to 
5 mm, but oxygen transfer efficiency decreases rapidly 
within this range (Figure 8). 

5. Comparison of Linear and Circular Bubble 
Plume Models 

[ 28 ) The primary difference between the linear and cir¬ 
cular bubble plume models is the plume geometry. For a 
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Plume Water Flow Rate (m 3 /s) Bubble Radius (m) 


Figure 6, In-plume profiles predicted by linear bubble plume model, represented as (a-c) solid lines 
and (d-f) solid lines and symbols. Input data collected from Spring Hollow Reservoir, Virginia during 
diffuser operation with compressed air (2 July 2003, crosses) and pure oxygen (17 August 2003, solid 
circles, and 23 October 2004, open circles). Measured average in-plume temperature, DO, and plume 
water density represented as symbols (Figures 6a-6c). 


given plume cross-sectional area, the perimeter of a linear 
plume is much greater than for a circular or round plume. 
The initial estimated plume area and perimeter of the linear 
diffuser in SHR for July 2003 was about 50 m 2 and 600 m, 
respectively (Table 3). This corresponds to an equivalent 
radius and perimeter of approximately 4 m and 25 m, 
respectively, for a circular plume. In this case, the perimeter 
of the linear plume is 24 times greater than for the circular 
plume. For both the linear and circular plume models, 
entrainment of ambient water is directly proportional to 
local plume perimeter, local plume water velocity, and a 
(Table 1). The larger perimeter of the linear plume greatly 
increases ambient entrainment, which contributes to the 
negative buoyancy of the plume and causes the plume to 
decelerate more rapidly. This results in a lower plume rise 
height compared with the circular plume. 

[ 29 ) The sensitivity of the linear bubble plume model to 
various parameters was comparable to results for the circu¬ 


lar bubble plume model analysis. Gas flow rate had the 
greatest effect on DMPR predictions by both models, and 
the linear plume model was less sensitive to initial bubble 
radius than the circular plume model was [Wuest et al, 
1992]. The latter may be due to differences between 
temperature boundary conditions, which caused plume 
dynamics in SHR to be more influenced by thermocline 
depth. Initial bubble size greatly affected oxygen transfer 
efficiency for both diffuser geometries, decreasing rapidly 
as the radii increased beyond 1 mm and 3 mm for the linear 
and circular plumes, respectively. Both the linear and 
circular plume model predictions for DMPR were relatively 
insensitive to Fr a and initial plume area (Figure 7 and Wilest 
et al. [1992]). However, circular model. predictions for 
DMPR were more sensitive to a. Overall, the linear model 
was less sensitive to input parameters than the circular 
model, but this insensitivity is probably due to the relatively 
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Figure 7. Effect of linear bubble plume model parameters on depth of maximum plume rise (DMPR). 
Standard values for each parameter are indicated by the vertical dashed lines. Input data collected from 
Spring Hollow Reservoir, Virginia, during diffuser operation with compressed air (2 July 2003, crosses) 
and pure oxygen (17 August 2003, solid circles, and 23 October 2004, open circles). 


homogeneous boundary conditions in the hypolimnion of 
SHR caused by diffuser mixing. 

6. Summary and Conclusions 

[ 30 ] In the current work, the linear bubble plume model 
of McGinnis et al. [2001] was improved, and the updated 
model was evaluated using data collected from a full-scale 
hypolimnetic oxygenation system installed in Spring 
Hollow Reservoir, Virginia. Three diffuser experiments 
were conducted using compressed air and pure oxygen over 
a range of flow rates. Predicted plume rise height, spread¬ 
ing, and constituent profiles were compared with experi¬ 
mental observations. For July 2003 and October 2004, the 
DMPR was simulated well by the model. However, the 
model underestimated the plume rise height for August 

2003, when the gas flow rate was comparatively low. The 
model underpredicted the final plume temperature by ap¬ 
proximately 0.3° and 0.2°C for July 2003 and October 

2004, respectively. The model predicted the plume DO 



Figure 8. Effect of initial bubble radius on oxygen transfer 
efficiency predicted by linear bubble plume model. Input 
data collected from Spring Hollow Reservoir, Virginia, 
during diffuser operation with compressed air (2 July 2003) 
and pure oxygen (17 August 2003 and 23 October 2004). 
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profiles very well for ail three diffuser tests, including 
simulating the initial rapid transfer of oxygen immediately 
above the diffuser. 

[ 31 ] A sensitivity analysis was performed to determine 
the effect of various model input parameters. The DMPR 
and induced water flow rate were most influenced by gas 
flow rate and initial bubble radius, moderately sensitive to a 
and the initial plume area, and insensitive to Fr a , A, and 
ambient dissolved nitrogen. Oxygen transfer within the 
hypolimnion was independent of all parameters except 
initial bubble radius, decreasing from nearly 100% to 
around 20% for air and pure oxygen as the initial radii 
increased from approximately 1 mm to 1 cm. 

[ 32 ] The linear bubble plume model for hypolimnetic 
oxygenation has been successfully validated. The model 
can be used to design lake and reservoir oxygenation 
systems and to optimize existing systems to maximize 
oxygen addition. Application of the linear plume model 
requires knowledge of the gas flow rate, initial bubble 
radius, initial plume area, and near-field constituent profiles 
(boundary conditions). Additionally, several empirical 
parameters must be estimated including a and A. Even 
though all of the model inputs may not be known with 
certainty for a given aeration or oxygenation system, the 
model can be used for preliminary design and coarse 
optimization. Also, plume rise height, water flow rate, and 
oxygen transfer efficiency were found to be primarily 
dependent on gas flow rate and initial bubble radius. 

[ 33 ] Operation of bubble plumes for hypolimnetic oxy¬ 
genation usually alters the ambient temperature and DO 
conditions of a water body. Plume dynamics and oxygen 
transfer are strongly related to the near-field water column 
properties, establishing a feedback loop that continually 
changes plume dynamics. This complex plume-lake inter¬ 
action should be accounted for in the design and operation 
of bubble plume diffusers. McGinnis et al. [2001] per¬ 
formed a preliminary coupling of the linear bubble plume 
model with an existing reservoir model, CE-QUAL-W2 
[Cole and Wells , 2003], and obtained encouraging results. 
Efforts are currently under way to further develop the 
coupled model to predict plume performance and the near- 
and far-field reservoir responses. In the absence of near- 
field boundary profiles, far-field or simulated constituent 
profiles can be used to provide a reasonable estimate of 
plume performance. 

Notation 

b plume radius (circular bubble plume), m. 

d bubble diameter, mm. 

C dissolved concentration, mol/m 3 . 

E entrainment factor, nr/s. 

F d dissolved species flux, mol/s. 

F a gaseous species flux, mol/s. 

F s salinity flux, kg/s. 

F t temperature flux, °C m 3 /s. 

Fr Froude number [ ]. 

g gravitational acceleration, m/s 2 . 

H Henry’s constant, mol/m 3 /bar. 

K l mass transfer coefficient, m/s. 

L plume length, m. 

M water momentum, m 4 /s 2 . 


N number flux of bubbles, 1/s. 

P pressure, bar. 

Q plume flow rate, m 3 /s. 
q actual gas flow rate per unit diffuser length, 
m 2 /h. 

Ri Richardson number [ ]. 

R bubble radius, m. 

5 salinity, g/kg. 

T temperature, °C. 
v velocity, m/s. 

W plume width, m. 
y gaseous concentration, mol/m 3 . 

Z depth, m. 

Greek letters 

a entrainment coefficient [ ]. 

A spreading coefficient [ ]. 
p density, kg/m 3 . 

Subscripts 

3,2 Sauter-mean. 

G Gaussian profile. 

O oxygen. 

N nitrogen. 

T top-hat profile. 
a ambient water. 
b bubble. 

i gas species, oxygen or nitrogen. 
o initial. 

p plume water and gas mixture. 
w plume water. 
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ABSTRACT 

Weiss, R.F., 1974. Carbon dioxide in water and seawater: the solubility of a non-ideal gas. 
Mar. Chem., 2: 203-215. 

New measurements of the solubility of carbon dioxide in water and seawater confirm 
the accuracy of the measurements of Murray and Riley, as opposed to those of Li and Tsui. 
Corrections for non-ideal behavior in the gas phase and for dissociation in distilled water 
are required to calculate solubility coefficients from these sets of data. Equations for the 
solubilities of real gases are presented and discussed. Solubility coefficients for carbon 
dioxide in water and seawater are calculated for the data of Murray and Riley, and are 
fitted to equations in temperature and salinity of the form used previously to fit the solu¬ 
bilities of other gases. 


INTRODUCTION 

Until very recently, direct determinations of the solubility of C0 2 in sea¬ 
water have been limited to a few measurements of Krogh (1904). Most authors 
have preferred to use approximations based on the solubility of C0 2 in aqueous 
sodium chloride solutions. Buch et al. (1932) prepared tables of seawater C0 2 
solubilities from the Bohr (1899) data for sodium chloride solutions, using 
the assumption that the effect of sea salt on the solubility is equal to the ef¬ 
fect of an equal weight of sodium chloride. Lyman (1957) refined this assump¬ 
tion by considering the salting-out of the various constituents of sea salt, sug¬ 
gesting that the Buch et al. values may be about 0.5% too high. 

The lack of direct measurements of C0 2 solubility in seawater has recently 
been alleviated by two independent sets of data: those of Li and Tsui (1971), 
determined by infrared analysis, and those of Murray and Riley (1971), de¬ 
termined gravimetrically. Unfortunately, the agreement between these sets of 
measurements is poor. Whereas the data of Li and Tsui support the values as¬ 
sumed by Buch et al., the solubility data of Murray and Riley are as much as 
3.8% lower than those of Li and Tsui at higher temperatures and salinities. 

The agreement is better at lower temperatures and in distilled water, so that 
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the differences cannot be expressed either as a constant offset or as a con¬ 
stant factor. 

In hopes of resolving these discrepancies, it was decided to make several 
independent measurements of C0 2 solubility. If it could be shown that a few 
accurate measurements covering a wide range of temperature and salinity 
were consistently in close agreement with one of the two sets of data, this 
would provide strong support for the accuracy of that set at all measured 
temperatures and salinities. 

The second objective of this work was the fitting of the best C0 2 solubility 
data to an equation in temperature and salinity of the form used previously 
to fit the solubilities of N 2 , 0 2 , and Ar (Weiss, 1970), He and Ne (Weiss, 
1971a), H 2 (Crozier and Yamamoto, 1974), and Kr (Weiss and Kyser, in prep¬ 
aration). C0 2 is many times more soluble than any of the gases treated previ¬ 
ously, and in the gas phase its departures from the ideal gas approximation 
are large compared to the accuracy with which its solubility can be measured. 
Therefore, before the use of such equations is justified, it is necessary to use 
Henry’s law for real gases and to test its validity over a range of partial pres¬ 
sures of C0 2 . Also, it is necessary to test the validity of the logarithmic Set- 
ch6now salting-out relation (Weiss, 1970,1971b) for C0 2 in seawater. 

Both Li and Tsui (1971), and Murray and Riley (1971), noted that their 
measurements would be affected by the dissociation of dissolved hydrated 
C0 2 to form bicarbonate, and accordingly they acidified their seawater sam¬ 
ples to a sufficiently low pH so that this effect could be neglected. However, 
both studies failed to take this effect into account for their distilled-water 
measurements, which were made without acidification. Although this effect 
cancels in direct comparisons between the two sets of measurements, the 
proper use of Henry’s law and of the Setch^now relation requires that the 
data be corrected for dissociation. 

All distilled water C0 2 solubility data discussed in the following sections 
have been corrected for dissociation using the following approximation: 

[co 2 ] = [sco 2 ] - Vk.'IscOj~] (1) 

where [C0 2 ] is the sum of the concentrations of. dissolved C0 2 and undis¬ 
sociated hydrated C0 2 , [2C0 2 ] is the sum of the concentrations of all dis¬ 
solved C0 2 species as measured in the solubility experiment, and K{ is the 
first apparent dissociation constant at the temperature of the measurement.* 1 
Solubilities measured in distilled water for Pco, — a ^ m a *’ e thus reduced 
by - 0.18% at 0°C, to ~ 0.46% at 40°C. 

HENRY'S LAW AND REAL GASES 

King (1969, ch.4) presents an excellent and thorough discussion.of Henry’s 
law and its application to real gases over a wide range of pressures. By ex- 


*’ Values of K\ were calculated according to Harned and Davis (1943, eq.15): 
log, 0 Kj = -3404.71/T + 14.8435-0.032786 T. 
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pressing the activity of the solute in the gas phase by its fugacity, and by 
making the appropriate thermodynamic correction for the expansion of the 
solution by the dissolved gas, the “modified Henry’s law equation” is obtained: 

fi = QtXi exp(Py,/ET) W 

where f { is the fugacity of gas i, x f is the mole fraction of i in solution, v t is 
the partial molal volume of i in solution, R is the gas constant, T is the abso¬ 
lute temperature, P is the total pressure at the liquid-gas-interface, and Q h 
the modified Henry’s law constant for i, is a function only of the temperature 
and the nature of the solvent. 

Conveniently, King uses the C0 2 -water system as an example to illustrate 
the application of eq.2. Using literature data, with special emphasis on the 
work of Wiebe and Gaddy (1939; 1940), Kang (1969, pp.219 220) demon¬ 
strates that eq.2 holds over the entire C0 2 -pressure range from zero to several 
hundred atmospheres, at a number of different temperatures.* 2 The use of 
the modified form of Henry’s law to calculate C0 2 solubilities over wide 
ranges of total pressures and C0 2 partial pressures, based on accurate solubil¬ 
ity measurements made at pressures near 1 atm, is therefore well justified.* 3 
This is particularly germane to studies dealing with the solubility of natural 
levels of atmospheric C0 2 . 

In treating the solution of atmospheric gases in natural waters, it is con¬ 
venient to express the solubility in terms of the Bunsen coefficient 0, thereby 
avoiding the problem of evaluating the mole fraction of the gas in mixed sol¬ 
vents such as seawater. The Bunsen coefficient is defined here as the volume 
of gas (STP) absorbed per unit volume of the solution, at the temperature of 
the measurement, when the total pressure and the fugacity are both 1 atm. 

This is preferable to the usual definition of 0 for ideal gases, in which the 
partial pressure is set at 1 atm and the total pressure is not specified, because 
the total pressure is required to define the system (see eq.2). Thus, C,-, the 
volume (STP) of gas i dissolved in a unit volume of solution at the tempera¬ 
ture of the measurement, is given by the relation: 

C,- = 0, fi exp[(l - P)vilRT] ( 3 > 

where: 

0i - Q: 1 vf pN exp(—i?,-/RT) ( 4 ) 

and vj is the volume of one mole of the real gas i at STP, p is the density of 
the solution, N is the number of moles in a unit weight of solvent, and P is in 
atmospheres. Since the solubility of most gases {v at 30 em 3 /mole) is decreased 


* 3 Note that helium solubility measurements by these workers are in excellent agreement 
with recent microgasometric data (Weiss, 1971a). 

* 5 A more detailed discussion of CO, solubility as a function of pressure is given m the 

** Pressures are given here in atmospheres because of the widespread use of this unit to 
define standard states. Atmospheres may be converted to bars (10 5 newtons/m ), the more 
fundamental unit commonly used in high-pressure work, by multiplying by 0.986923. 
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by only ~0.15% per atmosphere increase in total pressure, the small deviations 
from 1 atm total pressure which are encountered in many natural conditions 
may be neglected, thereby reducing the exponential term in eq.3 to unity. 

In the special case of C0 2 , solubilities are generally used in calculations of 
chemical equilibria and are therefore best given in molar (moles/1 of solution) 
or gravimetric (moles/kg of solution) units: 

[C0 2 ] = K 0 f COi exp[(l - P)U C0 JRT] (5) 

where the constant K 0 equals p/V* in molar units (moles/1 • atm) or (3/p in 
gravimetric units (moles/kg * atm). Again, the exponential term may be taken 
as unity at total pressures near 1 atm. 

In the following discussions, gas volumes and fugacity-pressure corrections 
are based on the virial equation of state. The appropriate expansions arc carried 
only as far as the second virial coefficient B(T), since the termscontaining 
the third virial coefficient were always found to be negligible (~10' s ). Values 
of B(T) in cm 3 /mole, for C0 2 (Sengers et al., 1971) in the range 265-320°K, 
are well represented by a power series: 

B(T) =-1636.75 + 12.0408 T~ 3.27957 • 10" 2 T 2 + 3.16528 • 1(T S T 3 (6) 

Thus, for the pure gas (neglecting the small contribution of water vapor), the 
solubility measurements discussed in the following sections have been corrected 
using the approximations (Guggenheim, 1967, pp:91 and 97): 

V{P,T) = V* (P,T) + B(T) ( 7 ) 

and: 

\n(f/P) = B(T)P/RT W 

where V is the volume of one mole of the real gas, and V* is the volume of 
one mole of ideal gas. 

In most natural applications which do not require accuracies greater than 
-0.7%, the fugaeity in eq.4 may be taken as equal to the partial pressure. 
However, if greater accuracy is desired, the fugaeity must be calculated. For 
a nearly pure gas phase, eq.8 will suffice, but for multi-component gas phases, 
such as C0 2 in air, it is necessary to calculate fugaeity in the mixture. 

Useful equations for the calculation of fugacities in binary mixtures are 
given by Guggenheim (1967, pp.175—177): 

f-y = XyP exp[(J3jj + 2x 15 12 ) F/RT] (9) 

where the subscripts 1 and 2 refer to the two components of the mixture, x is 
the mole fraction, and B x \ is the second virial coefficient of the pure gas 1. 

The quantity 6 12 is defined by: 

^12 = 2 + ®22) + 5 12 


( 10 ) 
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where B 12 is the cross virial coefficient for interactions between gas 1 and gas 
2 molecules, and B n and B 22 refer to the pure gases 1 and ^respectively. In 
addition to the binary mixtures, equations 8 and 9 may also be applied to 

mixtures of C0 2 in air when x COj « 1. 

Few direct measurements of cross-vinal coefficients exist, so that the quan¬ 
tity 5 jo must be evaluated from theoretical considerations. The quantity 
2>CO _ a j r has been evaluated from the Lennard-Jones (6-12) potential model 

following the method of Hirschfelder, Curtiss and Bird (1954,^.166-170). 
Hie results in cm 3 /mole for the temperature range 273 to 313°K are well 
represented by the linear relation: 

5 C0 2 -air= 57.7-0.118 T < U) 

The magnitude of this correction, as opposed to assuming 6=0 (Lewis and 
Randall rule), is about 0.2% in the overall solubility calculation, or about one 
third of the total deviation from ideal gas behavior. 

Calculations of C0 2 fugacity using the virial equation of state are valid to 
within 0.1% for total pressures up to ~10 atm. Higher pressures require the 
use of more sophisticated equations of state, such as the Benedict, Webb and 
Rubin (1940; 1942) equation, which may be used to at least 500 atm (see 
Appendix). 

DATA FITTING 

The corrected data of Murray and Riley (1971) were fitted to the same 
equation in temperature and salinity which has been used to fit the solubili¬ 
ties of many other gases (see Introduction). This equation is derived from the 
integrated van’t Hoff equation and the logarithmic Setch&iow salinity depen¬ 
dence, and has the form (Weiss, 1970): 

In K 0 = A x + A 2 (100/T) + A 3 In(T/100) + 

<S%)[B 1 + B a (T/100) + BAT/100) 2 ] (12) 

where K<> may be expressed either in moles/1 • atm, referring to a liter of solu¬ 
tion at the temperature of the measurement and an atmosphere fugacity in 
the gas phase, or in moles/kg • atm, referring to a kilogram of solution. The 
A’s and B's are constants, T is the absolute temperature, and is the salinity 
in parts per thousand. The method of fitting the data, and the conditions for 
justification of the Setch6now equation and the number of terms used in the 
expansion, are discussed elsewhere (Weiss, 1970). 

The corrected Murray and Riley C0 2 data showed good agreement with 
the equations at all stages of the fitting procedure. The data show a root-mean- 
square deviation from the final fitted equation of 1.4 • 10 4 moles/1 • atm in 
K 0 , or about 0.3%. Agreement with the Setch^now relation is shown in Fig.l 
by the random nature of the deviations as a function of salinity for each mea¬ 
sured isotherm. The fitted values of the constants in eq.12 are listed in Table I 
for K 0 in molar and in gravimetric units. Values of K 0 at various temperatures 
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Fig.l. Deviations of the Murray and Riley (1971) CO, solubility data (corrected for dis¬ 
sociation and non-ideality) in moles/1 • atm from the fitted curve (eq.l 2), plotted against 
salinity. Measurements’ made at constant temperature are connected by lines, and are 
labeled to the nearest °C. The dashed lines show ±1 standard deviation. 


aad salinities calculated from these constants are listed in Tables II and III. 

Salting-out constants (the J3’s in eq.12) were also determined for the solubil 
ity of CO i in aqueous NaCl solutions using the data of Bohr (1899), corrected 
for dissociation. With NaCl concentration expressed as weight percent in solu- 

TABLEi 

Constants for the calculation of the solubility of CO, in molar and gravimetric units 

according to eq.l2 


Units of K 0 

mol w'I ■ atm moles/kg • atm 


A . 

-58.0931 

-60.2409 

a 7 

90.5069 

93.4517 

A 3 

22.2940 

23.3585 

3, 

0.027766 

0.023517 


-0.026888 

-0.023656 


0.0050678 

0.0047036 
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TABLE II 

The solubility coefficient K 0 (10' J moles/l • atm) 


Salinity (°4>) 

T(°C) 0 10 20 30 34 35 36 


-1 

— 

0 

7.758 

1 

7.468 

2 

7.174 

3 

6.905 

4 

6.650 

5 

6.408 

6 

8.178 

8 

5.751 

10 

5.366 

12 

5.017 

14 

4.700 

16 

4.412 

18 

4.149 

20 

3.910 

22 

3.691 

24 

3.491 

26 

3.307 

28 

3.138 

30 

2.983 

32 

2.840 

34 

2.708 

36 

2.587 

38 

2.474 

40 

2.370 


— 

7.273 

7.364 

6.990 

7.081 

6.723 

6.813 

6.469 

6.558 

6.229 

6.317 

6.001 

6.088 

5.785 

5.871 

5.580 

5.469 

5.200 

5.105 

4.867 

4.776 

4.646 

4.477 

4.264 

4.205 

4.008 

3.958 

3.775 

3.732 

3.562 

3.526 

3.368 

3.337 

3.190 

3.164 

3.027 

3.005 

2.878 

2.859 

2.741 

2.725 

2.615 

2.601 

2.498 

2.487 

2.391 

2.382 

2.292 

2.284 

2.201 


6.903 

6.760 

6.635 

6.498 

6.382 

6.261 

6.143 

6.017 

5.916 

5.795 

5.701 

5.585 

5.497 

5.386 

5.303 

6.196 

4.945 

4.846 

4.621 

4.529 

4.327 

4.243 

4.062 

3.983 

3.820 

3.747 

3.600 

3.533 

3.400 

3.337 

3.217 

3.158 

3.050 

2.995 

2.897 

2.846 

2.756 

2.709 

2.627 

2.583 

2.509 

2.468 

2.400 

2.361 

2.299 

2.263 

2.207 

2.173 

2.121 

2.090 


6.724 

6.689 

6.465 

6.431 

6.219 

6.187 

5.986 

5.955 

5.766 

6.736 

5.557 

5.528 

5.358 

5.331 

5.170 

5.144 

4.822 

4.797 

4.507 

4.485 

4.222 

4.201 

3.964 

3.945 

3.729 

3.712 

3.516 

3.499 

3.322 

3.306 

3.144 

3.130 

2.982 

2.968 

2.833 

2.821 

2.697 

2.685 

2.572 

2.561 

2.457 

2.447 

2.352 

2.342 

2.254 

2.246 

2.165 

2.157 

2.082 

2.074 


38 

40 

6.620 

6.551 

6.364 

6.298 

6.123 

6.060 

5.894 

5.833 

5.677 

5.619 

5.472 

5.416 

5.277 

5.223 

5.092 

5.040 

4.749 

4.702 

4.440 

4.396 

4.160 

4.119 

3.906 

3.869 

3.676 

3.641 

3.466 

3.434 

3.275 

3.245 

3.101 

3.073 

2.942 

2.915 

2.796 

2.771 

2.662 

2.639 

2.540 

2.518 

2.427 

2.407 

2.323 

2.305 

2.228 

2.211 

2.140 

2.124 

2.059 

2.044 


tion, and K 0 in moles/1 • atm, the constants are: B t = -0.68330, S 2 - 0.40911, 
B 3 = -0.064989. Averaged over the temperature range 0—40 C, K 0 is 0.8% 
higher in a 3.6% Nad solution than K 0 in seawater of 36%o salinity, calculated 
from the salting-out constants in Table I. This difference probably lies within 
the error of the Lyman (1957) prediction that K 0 in the NaCl solution would 
be 0.5% higher. 


EXPERIMENTAL METHOD 

Measurements of CO* solubility were carried out by the microgasometric 
technique used previously to measure He and Ne solubilities (Weiss, 1971a). 
High-purity CO* (certified > 99.99%) was supplied by Matheson Gas Products 
andgas chromatographic analysis showed < 0.01% air contamination. Because 
of the high solubility of C0 2 compared to He and Ne, the amount of degassed 
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water added to the equilibration chamber was reduced to about one third of 
the total gas volume. Under these conditions, equilibration was extremely 
rapid (first-order rate constant r s 30 sec) so that the 6 10 minutes allowed 
for equilibration were more than adequate. Permeation of C0 2 through the 
indicator drop was rapid, and required that the chamber above the drop be 
filled with pure C0 2 prior to each equilibration. With this procedure, no drift 
in the indicator drop could be detected over a period of 1 h. 

In order to obtain the most direct comparisons with the work of Murray 
and Riley and of Li and Tsui, chemical procedures similar to their’s were fol¬ 
lowed. Distilled water measurements were performed without acidification 
and the corrections for dissociation were made. Surface seawater, collected 
at La Jolla and evaporated to increase its salinity by ~2%o, was passed through 
a 0.45 m filter, poisoned with 1 mg/liter HgClj, and its salinity determined to 
+ 0.004°/oo by an inductive salinometer. Sulfuric acid (~2 N) was then added 
to bring the pH to 2.2 and the salinity was adjusted by gravimetric determina¬ 
tion of the amount of H 2 0 added in the acid solution. Following Murray and 


TABLE III 


The solubility coefficient K 0 (10 "* moles/kg * atm) 


Salinity (7«) 


T(°C) 

0 

10 

20 

30 

—1 

— 

— 

7.158 

6.739 

0 

7.758 

7.305 

6.880 

6.479 

1 

7.458 

7.024 

6.616 

6.232 

2 

7.174 

6.758 

6.367 

5.999 

3 

6.904 

6.506 

6.131 

5.777 

4 

6.649 

6.267 

5.907 

5.568 

5 

6.407 

6.040 

5.695 

5.369 

6 

6.177 

5.825 

5.493 

5.180 

8 

5.752 

5.427 

5.120 

4.831 

10 

5.367 

5.067 

4.784 

4.516 

12 

5.019 

4.741 

4.479 

4.231 

14 

4.703 

4.446 

4.202 

3.972 

16 

4.416 

4.177 

3.951 

3.738 

18 

4.155 

3.933 

3.723 

3.524 

20 

3.916 

3.710 

3.515 

3.330 

22 

3.699 

3.507 

3.325 

3.152 

24 

3.499 

3.321 

3.151 

2.990 

26 

3.317 

3.150 

2.992 

2.841 

28 

3.149 

2.994 

2.846 

2.705 

30 

2.995 

2.850 

2.712 

2.580 

32 

2.854 

2.718 

2.589 

2.466 

34 

2.723 

2.596 

2.476 

2.360 

36 

2.603 

2.484 

2.371 

2.263 

38 

2.492 

2.381 

2.275 

2.174 

40 

2.389 

2.285 

2.186 

2.091 


34 

35 

36 

38 

40 

6.579 

6.539 

6.500 

6.422 

6.345 

6.325 

6.287 

6.249 

6.175 

6.101 

6.085 

6.048 

6.012 

5.941 

5.870 

5.857 

5.822 

5.788 

5.719 

5.651 

5.642 

5.608 

5.575 

5.509 

5.444 

5.438 

5.40.5 

5.374 

5.310 

5.248 

5.244 

5.213 

5.182 

5.122 

5.062 

5.060 

5.031 

5.001 

4.943 

4.885 

4.720 

4.693 

4.666 

4.612 

4.658 

4.413 

4.388 

4.363 

4.313 

4.263 

4.136 

4.112 

4.089 

4.042 

3.997 

3.884 

3.862 

3.840 

3.797 

3,755 

3.655 

3.635 

3.615 

3.575 

3.536 

3.448 

3.429 

3.410 

3.373 

3.336 

3.258 

3.241 

3.223 

3.189 

3.154 

3.086 

3.069 

3.053 

3.021 

2.989 

2.928 

2.912 

2.897 

2.867 

2.837 

2.783 

2.769 

2.755 

2.727 

2.699 

2.651 

2.638 

2.624 

2.598 

2.572 

2.530 

2.517 

2.505 

2.480 

2.455 

2.418 

2.406 

2.395 

2.372 

2,349 

2.316 

2.305 

2.294 

2.272 

2.250 

2.221 

2.211 

2.201 

2.180 

2.160 

2.134 

2.125 

2.115 

2.096 

2.077 

2.054 

2.045 

2.036 

2.018 

2.000 
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Riley (1971), the effect of the added H 2 S0 4 molecules on the salinity was 
neglected as being insignificant with respect to altering the solubility of C0 2 . 
The water was degassed using the same vacuum extraction method as in the 
previous work (Weiss, 1971a). 

The results were calculated in the same manner as previous microgasometnc 
measurements, except that the corrections for non-ideal behavior discussed in 
the foregoing section were applied. The correction for the expansion of the 
aqueous phase during equilibration was made using a partial molal volume of 
32.3 cm 3 per mole of C0 2 (see Appendix). The overall accuracy of these C0 2 
solubility measurements is estimated as ±0.2%. 


RESULTS AND CONCLUSIONS 

The microgasometnc experimental values for K 0 in moles/1 • atm at several 
temperatures and salinities, corrected for dissociation and non-ideality, are 
given in Table IV. Deviations of the microgasometnc results from the fitted 
Murray and Riley data, plotted in Fig.2, show excellent agreement without 
systematic deviations. At each of the three different temperature and salinity 
conditions, the range of microgasometnc points brackets the fitted curve. The 
average deviation of all the measured points is 1.5 • 10 moles/1 • atm, and the 
root-mean-square deviation is 1.2’ 10" 4 moles/1’ atm. for the deviations 

of the data of Li and Tsui, taken from their equations (1) and (2a) after cor¬ 
rection for dissociation and non-ideality, are also plotted m Fig.2. Although 
their fresh-water measurements are in reasonably good agreement, the seawater 
data of Li and Tsui show large deviations from the fitted Murray and Riley 


TABLE IV 

Experimental C0 2 solubilities: microgasometnc determinations of K„ in moles/1 • atm 


Salinity (%o) 

t(°C) 

Ko’lO’ 

0.0 

20.61 

3.848 

0.0 

20.63 

3.843 

0.0 

20.60 

3.867 

0.0 

20.59 

3.834 

0.0 

20.61 

3.840 

35.330 

6.59 

5.071 

35.330 

6.60 

5.063 

35.330 

6.60 

5.068 

35.330 

6.69 

5.041 

35.330 

20.62 

3.245 

35.330 

20.63 

3.272 

35.330 

20.63 

3.245 

35.330 

20.63 

3.255 

35.330 

20.63 

3.264 

35.330 

20.63 

3.252 



212 



Pic 2 Deviations of various CO, distilled water and seawater solubility data in moles/l- atm 
fromtS fit to the corrected data (see text) of Murray and Riley (1971), plotted1 ag«nst 
temperature. The curves for Li and Tsui (1971) are taken from their equations (IJand^ 

(2a Wad like the data of Krogh (1904), have been corrected for non-ideality and dissocia 
tion’ the dashed portions of these curves represent the extrapolation of their equations 
beyond the range of their measurements. The magnitude of the error caused by 
for measurements made in distilled water at pCO, = 1 atm is shown by the Ime labeled HCO 


data. Fig.2 also shows the corrected data of Krogh (1904), which agree well 

with Murray and Riley at the higher temperatures. 

The microgasometric results provide independent confirmation of the accu¬ 
racy of the C0 2 solubility measurements of Murray and Riley (1971) and 
show the measurements of Li and Tsui (1971) to be in error by as.much_as 
~4% at the higher temperatures and salinities. The data of Murray and Ruey, 
corrected for the effects of non-ideal gas behavior and for dissociation in the 
f jj p t H lwl water measurements, are well represented by the Setch6now and 
integrated van’t Hoff equations used previously to fit the solubilities of several 
other comparatively ideal gases. Thus, C0 2 solubilities calculated from eq. , 
using the constants in Table I, are believed to provide the most accurate values 
in the literature, with an overall accuracy estimated at ±1 • 10 moles/1- atm 
or about 0.2%. The solubility of C0 2 obeys the modified form of Henry s law 
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(eq.2) and it is therefore necessary to take account of the total pressure, as 
well as to calculate the fugacity of C0 2 in the gas phase, if full use is to be 
made of the accuracy of these solubility data. 
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APPENDIX 

Carbon dioxide solubility at high pressures 

At pressures greater than ~ 10 atm, CO, fugacities calculated from the virial equation 
of state, even when carried to the third virial coefficient, show large deviations from the 
literature values (Deming and Deming, 1939). Far better results are obtained with the 
Benedict, Webb and Rubin equation of state for pure substances (Benedict et al., 1940) 
and for mixtures (Benedict et al., 1942), which fit the observed p—V—T data up to densi¬ 
ties of twice the critical density. In the following discussion, Benedict, Webb and Rubin 
equation constants for pure CO, are taken from the Bishnoi and Robinson (1971)* 5 
volume-dependent fit to the compressibility data of Reamer et al. (1944). Because the 
Benedict, Webb and Rubin equation is an explicit function of V and T, the solution for 
specific values of P is obtained by the method of successive approximations. 

Fig.3 shows the distilled water CO, solubility data of Wiebe and Gaddy (1940) for the 
temperature range 12—40°C and the pressure range 25—500 atm, plotted as In (f/x) against 
the total pressure P. Fugacities were calculated assuming a pure CO, phase at a pressure of 
P minus the vapor pressure of water. Negligible error was introduced by this approximation 
because of the small fraction of water vapor in the CO,-rich phase. The straight lines plotted 
at each of the six experimental temperatures were fitted to the Wiebe and Gaddy data by 
the method of least-squares. 

According to eq.2, if v is independent of P, then In (f/x) plotted as a function of P at 
constant T should give a straight line with intercept In Q and slope (D/RT). This relationship 
is strongly supported by the experimental data as shown in Fig.3. The root-mean-square 
deviation of the Wiebe and Gaddy measurements about the six straight lines is ~0.8%, and 
there is no systematic indication of departure from linearity. The slopes of these lines cor¬ 
respond to remarkably constant values of 5, with a mean of 32.3 cm 3 /mole and a standard 
deviation of 0.5 cm’/mole. Thus, the data show no significant variation in v, either as a func¬ 
tion of T over the range 12—40°C, or as a function of P up to 500 atm. 

Results of the fit to the corrected Murray and Riley distilled water data (eq.12), mea¬ 
sured at 1 atm, are also plotted in Fig.3. The overall agreement with the Wiebe and Gaddy 
data is good, although there may be a systematic difference at 12°C (the Wiebe and Gaddy 
data also show the greatest scatter at this temperature). On the average, the solubilities at 
1 atm obtained from the linear fits to the Wiebe and Gaddy data are 0.4% lower in (3 (or 
higher in Q), and show a 2% root-mean-square difference. 


* s The constants, converted from English to metric units (atm, 1, mole, K): A„ 1.9521940, 
B 0 = 0.033065171, C, = 170449.55, a « 0.25264437, b - 6.6041298-10~ 3 , c = 19692.012, 
a - 4.7117009-10"’, 7 = 4.3414415-10" 3 , R = 0.08205601. 
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Fig. 3 . High-pressure distilled water CQ, solubility data of Wiebe and Gaddy (1940) plotted 
as the logarithm of the fugacity to mole fraction ratio.,against the total pressure. Solubili¬ 
ties at 1 atm calculated from eq.12 are plotted for the temperatures of the Wiebe and 
Gaddy measurements. 


In summary, the high-pressure solubility data show that the modified form of Henry’s law 
(eqs.2, 3, or B) is valid for GO, in water over the entire pressure range of. 0-BOO atm, which 
corresponds to a dissolved CO, concentration range of 0-1.7 molar. Assuming that v for 
CO, is the same in seawater as in distilled water, as was shown by Ennset.ah (1966) for Q„ 
the solubility of CO, in seawater at high pressures may also be calculated from these equations. 
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The indices in the original Johnson and Braida correlation and the point of 
intersection between the two linear regions were adjusted to improve the agree¬ 
ment with all the data meeting the above criteria. The resulting correlation 

(G 12 ) is: 

J = 0.94H 0 - 757 (2 < H < 59.31 17-5) 
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Technology Center/Art Unit: 3752 
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For: METHOD AND APPARATUS FOR 
REDUCING THE INTENSITY OF 
HURRICANES AT SEA BY DEEP¬ 
WATER UPWELLING 


Customer No.: 20350 
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Commissioner for Patents 
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DECLARATION UNDER 37 C.F.R. § 1.132 


I, Isaac Ginis Ph.D., declare the following: 


1. I am a tenured Professor of Oceanography at the University of Rhode Island. I hold a 
doctorate in Geophysics, specializing in Oceanography and Meteorology. A copy of my resume 
is attached as Exhibit A. 

2. I am a leading expert in numerical modeling and forecasting of air-sea interaction during 
hurricanes. I have published over 70 papers in scientific journals and books on this topic and 
authored a chapter on hurricane-ocean interaction for the book “Global Perspectives on Tropical 
Cyclones” published by the World Meteorological Organization, Geneva, Switzerland in 1995. I 
have delivered over 100 lectures worldwide at national and international meetings as well as at 
oceanographic and atmospheric institutions and have advised 12 doctoral students and 
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postdoctoral scientists. The National Science Foundation, the National Oceanic and 
Atmospheric Administration, the U.S. Office of Naval Research supported my research. 

3. I have been leading the effort toward improvements of the GFDL/URI coupled hurricane 
model related to the ocean and wave coupling. This work involves close collaboration between 
my research group at URI and scientists at the NOAA’s National Centers for Environmental 
Prediction (NCEP) and Geophysical Fluid Dynamics Laboratory (GFDL). My research group 
has made a successful conversion of a research coupled hurricane-ocean research model to a 
fully automated real-time prediction system. This significant technological and computer 
programming effort has resulted in implementation of the GFDL/URI coupled hurricane model 
to operational forecasting at National Weather Service in 2001. This model is used by the 
National Hurricane Center for issuing official hurricane warnings. I continue to be responsible 
for maintaining and improving the coupled hurricane operational forecast system at NCEP. 

4. I have been retained by the assignee of the above-identified application to assist in 
responding to the Office Action. I have no financial interest in the assignee or the outcome of 
this patent application, including whether it issues as a patent or not. I am being compensated 
for the time spent on this matter at the rate of $200/hr, plus reasonable expenses. 

5. I have read and am familiar with the above-referenced application. I have also read and 
am familiar with the Office Action mailed May 4, 2007 ("Office Action") pertaining to this 
application. 

6. It is my understanding that the claims currently under examination, which relate to, inter 

alia, methodologies for reducing the intensity of a hurricane, were rejected as allegedly being 

wholly inoperative, lacking credible utility, and not enabled. Specifically, the Examiner asserts: 

.. .applicant admits in his arguments, "submersibles of the kind 
required for this application do not presently exist." It seems that 
applicant wishes that someone will come along and develop the 
technology required to make the required submersibles, thereby 
enabling the present invention. Therefore, it is impossible for one 
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of ordinary skill in the art at his time to make and or use this 
invention. Office Action at page 6. 

7. The disclosure contains sufficient evidence and reasoning that a person of skill in the art 
would appreciate that since a hurricane draws energy from the heat content of the upper ocean, it 
is generally accepted that a large area of cooled ocean surface can suppress hurricane intensity. 
For example, the specification states: 


[0005] Because tropical storms draw their energy from the heat content 
of the upper ocean, it is generally accepted that a large area of cooled 
ocean surface can suppress hurricane intensity. Numerical modeling 
studies at the Massachusetts Institute of Technology suggests that 
reduction of sea surface temperature by 2.5°C in the storm’s central core 
would eliminate the thermodynamic conditions that sustain hurricanes. 
Other numerical model studies by independent researchers corroborate 
these results. In addition, analyses of measurements from past hurricanes 
show a strong correlation between lack of hurricane intensification and 
conditions that favor cold-water upwelling by the storm’s own winds, 
such as a shallow thermocline or slow forward speed. Finally, there is 
clear evidence that hurricanes weaken (or do not intensify under 
otherwise favorable conditions) when a hurricane crosses the cold 
“wake” of a previous storm. 


I agree with this assertion in the specification, which is consistent with the knowledge of one 
skilled in the art at the time of filing the application. Indeed, as detailed below, I conducted 
numeric modeling that illustrates that reduction of the temperature in the storm’s central core 
would reduce the wind speed of the hurricane, and hence the intensity of the hurricane. 


8. Variability in hurricane intensity originates from two sources: internal variability and 
environmental interactions. An important aspect of environmental interaction is the coupling 
between the storm and the underlying ocean. It is well known that tropical cyclones are driven by 
turbulent heat fluxes from the ocean (Ooyama, 1969) and it is natural to suppose that variability 
in surface conditions and in the response of the upper ocean to the passing storm plays a role in 
controlling storm intensity. 
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9. The effect of air-sea interaction as a negative feedback on tropical cyclone development 
and intensity has been well established. It is known that strong surface winds in a tropical 
cyclone induce turbulent mixing in the upper ocean and entrainment of the underlying cold water 
into the ocean mixed layer, which cools and deepens (e.g., Price 1981, Bender et al. 1993, Ginis 
2002). Both observational and real case numerical studies (e.g.. Black, 1983; Bender and Ginis 
2000) showed that the SST anomalies induced by tropical cyclones can reach up to 5-6°C. 
Studies also showed that tropical cyclone intensity is more sensitive to the local SST changes 
under the hurricane core than to those beyond the core area (e.g., Emanuel, 1999; Shen et al., 
2000). Therefore, it can be expected that cooling of the ocean area underneath the hurricane core 
may reduce its intensity. 

10. Numerical modeling studies suggest that a reduction of sea surface temperatures by 2.5° 
C in the storm’s central core would eliminate the conditions that sustain hurricanes (Emanuel, 
1986,1988). Numerical simulations with coupled models (Khain and Ginis, 1991, Bender et al., 
1993) and the operational GFDL/URI coupled hurricane-ocean model forecasts (Bender et al, 
2007) indicate that most tropical cyclones experience a noticeable reduction of intensity owing to 
their coupling with the ocean, though the degree of reduction depends on many aspects of the 
upper ocean thermal structure and the speed of translation, size and initial intensity of the 
tropical cyclone. This coupling also depends on the surface heat and momentum exchange, 
which affects mean currents and turbulence in the ocean mixed layer and thereby influences the 
degree of mixing through the seasonal thermocline. 

11. For this declaration, I have conducted a set of numerical experiments to assess the impact 
on hurricane intensity of a region with the reduced upper ocean heat content (hereafter “cooled” 
region) placed in front of a moving hurricane. I used the NOAA operational GFDL/URI coupled 
hurricane-ocean model (Bender et al. 2007) with idealized, but realistic, oceanic and atmospheric 
conditions. The GFDL hurricane prediction system became operational in 1995 as the U.S. 
National Weather Service’s official hurricane model. Since that time, it has provided forecast 
guidance to forecasters at the NWS’s Tropical Prediction Center (TPC) and has been the most 
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reliable forecast model for track error during the past decade (Table 1, courtesy of James 
Franklin, TPC). Table 1 is attached as Exhibit B. 

12. Important upgrades to the physics and spatial and vertical resolutions have been 
introduced to the GFDL/URI model over the years. These upgrades led to significant increase in 
GFDL/URI model intensity forecast skill (Bender et al, 2007). The operational GFDL/URI 
model presently has about 9 km spatial resolution in the innermost movable mesh. For this 
study, I increased the finest resolution to about 4.5 km to improve further its intensity forecast 
skill. 

13. For these idealized experiments, I integrated the coupled model for 120 hours starting 
with a normal size initial hurricane vortex embedded in specified initial, horizontally uniform 
environmental conditions. A GATE (Global Atlantic Tropical Experiment) III condition in the 
tropics is used for the atmospheric environmental thermal profile, which has air temperature of 
27°C and relative humidity of 84% at the lowest model level (about 35 m). A steady easterly 
environmental wind of 5 m/s was used (the hurricane is moving north-west). In the control 
experiment, the ocean was initially horizontally uniform and motionless with the sea surface 
temperature of 28.5°C. The initial vertical temperature profile with a mixed layer depth of 30 m 
is typical for the northern Gulf of Mexico in September. 

14. Simulation of the SST response in the control experiment is shown in Figure 1, attached 
as Exhibit C. A typical SST cooling pattern is generated due to hurricane-ocean interaction. The 
largest cooling is on the right sight of the hurricane track, consistent with observations and other 
modeling studies (Ginis, 2002). 

15. In the first two sensitivity experiments, I introduced the cooled regions in front of the 
moving hurricane as shown in Figure 2, attached as Exhibit D. The size of each region is about 
400 km in the along-track direction and covered the entire computational domain in the cross¬ 
track direction. The temperature anomalies in the regions are 1°C and 2°C, correspondingly, 
which are evenly distributed over the depth of the mixed layer. 
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16. The main results are summarized in Figures 3 and 4 (Figure 3 is attached as Exhibit E 
and Figure 4 is attached as Exhibit F). The temperature anomalies (SSTA) underneath the 
hurricane core, defined as a circular area around the storm center with i?=l 00 km, are shown in 
Figure 3. The SSTA are greatly reduced when the hurricane crosses the cooled regions. 
Evolution of hurricane central pressure and maximum winds in the numerical experiments are 
shown in Figure 4. In both sensitivity experiments the hurricane intensity was reduced after the 
storm encounters the cooled regions. The maximum winds were reduced from about 145 kts to 
about 135 kts (6% reduction) in the 1°C swath experiment and to about 130 kts (10% reduction) 
in the 2°C swath experiment. 

17. I conducted an additional experiment in which the size of the 2°C cooled region was 
doubled along the track direction. As result, the hurricane intensity was further reduced from 
about 145 kts to about 120 kts (22% reduction). These sensitivity experiments clearly indicate 
that both the size and magnitude of the cooled area encountered by a moving hurricane make 
important impact on the hurricane intensity reduction. 

18. A list of references cited herein are attached as Exhibit G. 
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19. All statements made herein of my own knowledge are true and that all statements made 
on information and belief are believed to be true; and further, that the statements were made with 
the knowledge that willful false statements and the like so made are punishable by fine or 
imprisonment, or both, under section 1001 of Title 18 of the United States Code, and such willful 
false statements may jeopardize the validity of the application or any patents issuing thereon. 



Isaac Ginis, Ph.D. 


Date JloKov 
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Curriculum Vita 


ISAAC GINIS 

Professor of Oceanography 

Graduate School of Oceanography 

University of Rhode Island 

Narragansett, RI02882 

Tel: (401) 874-6484; Fax: (401) 874-6728 

e-mail: iginis@gso.uri.edu 

WWW: www.DQ.gso.uri.edu/Numerical/ig 


Professional Preparation: 


Kabardino-Balkarian State University, Nalchik, Russia Mathematics M.S. 
Institute of Experimental Meteorology, Obninsk, Russia Geophysics Ph.D. 


1977 

1986 


Appointments: 


2003-present 
2006 - Fall 
1998-2003 
1993-1998 

1990-1993 

1977-1989 


Professor of Oceanography, University of Rhode Island 
Visiting Research Scholar, Princeton University, Princeton, NJ 
Associate Professor of Oceanography, University of Rhode Island 
Associate Marine Research Scientist, Adjunct Professor of Oceanography 
University of Rhode Island 

Visiting Research Scientist, Geophysical Fluid Dynamics Laboratory/NOAA 
Princeton University, Princeton, NJ 
Research Scientist, 

Institute of Mechanics and Applied Mathematics, 

Kabardino-Balkarian State University, Nalchik, Russia 


Awards 

National Oceanic and Atmospheric Administration 2001 Outstanding Scientific Paper Award 
National Science Foundation Atmospheric Sciences Division 2001 Highlighted Research Project 
National Oceanic and Atmospheric Administration 2002 Environmental Hero Award 
National Oceanographic Partnership Program 2002 Excellence in Partnering Award 


Journal Publications (last 3 years) 

Bender, M.A., I. Ginis, R. Tuleya, B. Thomas, T. Marchok, 2007: The operational GFDL coupled 
hurricane-ocean prediction system and a summary of its performance. Mon. Wea. Rev. In press. 

Moon, I., I. Ginis, and T. Hara, B. Thomas, 2007: Physics-based parameterization of air-sea momentum 
flux at high wind speeds and its impact on hurricane intensity predictions. Mon. Wea. Rev. 135, 2869- 
2878. 

Moon, I., I. Ginis, and T. Hara, 2007: Impact of reduced drag coefficient on ocean wave modeling under 
hurricane conditions, Mon. Wea. Rev. In press. 

Yablonsky, R. M., I. Ginis, 2007: Improving the initialization of coupled hurricane-ocean models by 
assimilating mesoscale oceanic features. Mon. Wea. Rev. In press. 

Fan Y., I. Ginis, T. Hara, II Ju Moon, 2007: Energy and momentum budget across air-sea interface. Part I: 
Steady Uniform Wind, J. of Atmos. Sci. Submitted. 

Fan Y., I. Ginis, T. Hara, II Ju Moon, 2007: Energy and momentum budget across air-sea interface. Part II: 
Hurricane Wind, J. of Atmos. Sci. Submitted. 



Falkovich, A., and I. Ginis, 2005: Ocean data assimilation and initialization procedure for the Coupled 
GFDL/URI Hurricane Prediction System. J. Atmos. Oceanic Technol., 22, 1918-1932. 

Ginis, I., A.P. Khain, E. Morozovsky, 2004: Effects of large eddies on the structure of the marine 
boundary layer under strong wind conditions, J. Atmos. Sci., 61. 3049-3064. 

Moon, I.-J., I. Ginis, and T. Hara, 2004: Effect of surface waves on air-sea momentum exchange. Part II: 

Behavior of drag coefficient under tropical cyclones, J. Atmos. Sci., 61,2334- 2348. 

Moon I.J., I. Ginis, T. Hara, E. J. Walsh, and H. L. Tolman, 2003: Numerical modeling of sea surface 
directional wave spectra under hurricane wind forcing, J. Phys. Oceanogr., 33, 1680-1706. 

Frolov, S.A., Sutyrin, G.G., and I. Ginis, 2004: Asymmetry of the stabilized Gulf Stream system. J. 
Phys Oceanogr., 34, 1087-1102. 

Knutson, T.R., R.E. Tuleya, I. Ginis, 2004: Impact of climate change on hurricane intensity as simulated 
using regional high-resolution models. Book Chapter. Columbia University Press. 

Synergistic Activities 

I’ve been leading the effort toward improvements of the GFDL/URI coupled hurricane model related to 
the ocean and wave coupling. This work involves close collaboration between my research group at URI 
and scientists at the NOAA’s National Centers for Environmental Prediction (NCEP) and Geophysical 
Fluid Dynamics Laboratory (GFDL). My research group has made a successful conversion of a research 
coupled hurricane-ocean research model to a fully automated real-time prediction system. This significant 
technological and computer programming effort has resulted in implementation of the GFDL/URI 
coupled hurricane model to operational forecasting at National Weather Service in 2001. This model is 
used by the National Hurricane Center for issuing official hurricane warnings. We continue to be 
responsible for maintaining and improving the coupled hurricane operational forecast system at NCEP. A 
number of major improvements have been made to the operation model for the 2002-06 hurricane seasons 
which led to significant improvements of the GFDL model forecast skill. Through ongoing collaboration 
with GFDL, Navy’s Fleet Numerical Meteorology Center (FNMOC) and Joint Typhoon Warning Center 
(JTWC) we are transitioning the new version of the GFDL/URI hurricane model to operations at 
FNMOC for global tropical cyclone forecasting. Since 2006 my research group has been actively 
involved in developing the new generation Hurricane Weather, Research and Forecast (HWRF) model in 
collaboration with EMC/NCEP scientists. We focus our effort on developing and improving the ocean 
and wave components of the coupled HWRF system. We also collaborate with scientists at the Korean 
Oceanography and Development Institute in developing a coupled tropical cyclone-ocean model for the 
western Pacific. 

I have developed new courses at URI that cover various numerical methods applied for solving the 
fundamental equations governing atmospheric and oceanic motions, marine geophysics, and biophysics. 

Research Collaborators within Last 48 Months 

M. Bender, R. Tuleya, T. Marchok (GFDL/NOAA), A. Khain (HIUJ), G. Sutyrin, L. Rothstein, T. Hara 
(URI), Il-Ju Moon (JU, Korea), H. Tolman, N. Surgi, C. Lozano (NCEP/NOAA), E. Walsh (NOAA/ETL 
and NASA/GSFC), K. Emanuel (MIT), Sok Kuh Kang (KORD1). 

Graduate Students and Postdoctoral Fellows Advised/Sponsored 

Sergey Frolov (PhD), Evan Robertson (MS), Minoru Kadota (MS), Yalin Fan (PhD), Dr. Clark Rowley, 

Dr. Weixing Shen, Dr. Ray Richardson, Dr. Il-Ju Moon, Richard Yablonsky(PhD - present), Zhitao Yu 
(PhD- present), Erica Clay (PhD- present), Seunghoun Lee (PhD - present), Lou Licate (MS - present). 

Dr. Biju Thomas (present). Dr. Yalin Fan (present). 
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VERIFYING 

NUMBER 

GIM~ 

GFSI 

UKMI 

NGPI 

TIME 

OF CASES 

- 




00 

2252 

JWV 

8. 

8. 

8. 

12 

2128 

$8.8 , 

42.8 

44.9 

43.3 

24 

1952 


73.7 

77.7 

74.2 

36 

1741 


104.4 

108.8 

105.4 

48 

1511 

|2;8 

135.2 

136.7 

137.4 

72 

1159 


208.0 

195.0 

206.1 


TABLE 1 Average track errors in nautical miles [nmj for all forecasts run in the Atlantic 
between 1996 and 2005 for the GFDL, NCEP’s GFS, UKMET and Navy’s NOGAPS models. 
Results are for the time interpolated models. Since numerical weather prediction models are 
generally not available to the forecasters in time to make their forecasts, a simple technique 
exists to take the model forecasted position and intensity, and adjust the forecast to apply to the 
current synoptic time and initial conditions. This adjustment is usually 6 or 12 hours, depending 
on the availability of the last model guidance. These adjusted versions are known for historical 
reasons as interpolated models which are generally indicated by the letter “I ” at the end of the 
name (e.g., GFDIfor the GFDL interpolated model). 














Exhibit C 



34 


Idealised Hunicane Simulation (Control) 
SST anomaly at 120h forecast 


30 




Exhibit D 




Figure 2. SST anomalies at 120 h in two sensitivity experiments, where cooled regions of 1°C 
(left panel) and 2°C (right panel) temperature anomalies are placed in front of the hurricane. 
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Average SST decrease under ihe cydorte (averaged over 150Km radius) 
Anonisiy from bach ground SST 



Figure 3. SST anomalies (SSTA) underneath the hurricane core region, defined as a circular 
area around the storm center with R=100 km, as function of time. The rectangular box indicates 
the time when the hurricane center crossed the region with the reduced upper ocean heat 
content. 
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Figure 4. Evolution of central pressure and maximum winds in four numerical experiments: 1) 
control (blue line), 2) 1°C anomaly (left panel, red), 3) 2°C anomaly (middiet panel, red,) and 4) 
the same as 1) and 2) except the size of the cooled swath in the along track direction is twice 
larger. The green rectangular boxes indicate the time when the hurricane center crossed the 
regions with the reduced upper ocean heat content. 
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